function [result,bm,bs,t]  = iflt_fast(image, varargin)
%function [result,bm,bs,t]  = iflt_fast(image, neighbors, radius, bins, mode)
narginchk(1,5);
nargoutchk(1,4);
nvarargin = length(varargin);
% iflt_fast(image, neighbors, radius, bins, mode)
% iflt_fast(image, spoints, bins, mode)
if nvarargin < 1
    neighbors = 8;
else
    neighbors = varargin{1};
end
useSp = numel(neighbors) > 1;
if useSp
    spoints = neighbors;
    neighbors = size(spoints,1);
    if nvarargin < 2
        bins = 64;
    else
        bins = varargin{2};
    end
    if nvarargin < 3
        mode = 'nh';
    else
        mode = varargin{3};
    end
    if ~isnumeric(bins) || ~isnumeric(neighbors)
        error('iflt_fast:InvalidArg','bins, neighbors should be numeric');
    end
    if mod(neighbors, 2) ~= 0
        error('iflt_fast:InvalidArg','neighbors must be multiple times of 2');
    end
else
    if nvarargin < 2
        radius = 1;
    else
        radius = varargin{2};
    end
    if nvarargin < 3
        bins = 64;
    else
        bins = varargin{3};
    end
    if nvarargin < 4
        mode = 'nh';
    else
        mode = varargin{4};
    end
    if ~isnumeric(bins) || ~isnumeric(radius) || ~isnumeric(neighbors)
        error('iflt_fast:InvalidArg','bins, neighbors, radius should be numeric');
    end
    if mod(neighbors, 2) ~= 0
        error('iflt_fast:InvalidArg','neighbors must be multiple times of 2');
    end
    spoints = zeros(neighbors, 2);
    a = 2*pi/neighbors;
    spoints(:,1) = -radius*sin((1:neighbors)*a);
    spoints(:,2) = radius*cos((1:neighbors)*a);
end

% neighbors
% bins
% mode
% spoints
[ysize xsize] = size(image);

tic
% Determine the dimensions of the input image.
d_image = double(image);

miny=min(spoints(:,1));
maxy=max(spoints(:,1));
minx=min(spoints(:,2));
maxx=max(spoints(:,2));

% Block size, each LBP code is computed within a block of size bsizey*bsizex
bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1;
bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1;

% Coordinates of origin (0,0) in the block
origy=1-floor(min(miny,0));
origx=1-floor(min(minx,0));

% Minimum allowed size for the input image depends
% on the radius of the used LBP operator.
if(xsize < bsizex || ysize < bsizey)
  error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)');
end

% Calculate dx and dy;
dx = xsize - bsizex;
dy = ysize - bsizey;

% Fill the center pixel matrix C.
C = image(origy:origy+dy,origx:origx+dx);
d_C = double(C);

res = zeros(dy+1, dx+1, neighbors);
sum_r = zeros(dy+1, dx+1);

for i = 1:neighbors
  y = spoints(i,1)+origy;
  x = spoints(i,2)+origx;
  % Calculate floors, ceils and rounds for the x and y.
  fy = floor(y); cy = ceil(y); ry = round(y);
  fx = floor(x); cx = ceil(x); rx = round(x);
  % Check if interpolation is needed.
  if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)
    % Interpolation is not needed, use original datatypes
    N = d_image(ry:ry+dy,rx:rx+dx);
  else
    % Interpolation needed, use double type images 
    ty = y - fy;
    tx = x - fx;

    % Calculate the interpolation weights.
    w1 = (1 - tx) * (1 - ty);
    w2 =      tx  * (1 - ty);
    w3 = (1 - tx) *      ty ;
    w4 =      tx  *      ty ;
    % Compute interpolated pixel values
    N = w1*d_image(fy:fy+dy,fx:fx+dx) + w2*d_image(fy:fy+dy,cx:cx+dx) + ...
        w3*d_image(cy:cy+dy,fx:fx+dx) + w4*d_image(cy:cy+dy,cx:cx+dx);
  end
  D = d_C - N;
  sum_r = sum_r + D.^2;
  res(:,:,i) = D;
end

sum_r = sum_r .^ 0.5;
% res <= I_norm
for i = 1:neighbors
    res(:,:,i) = res(:,:,i) ./ (sum_r+eps);
end

I_norm_1 = res(:,:,1:2:end);
I_norm_2 = res(:,:,2:2:end);
y_high = (I_norm_1 + I_norm_2) / 2^0.5;
y_low = (I_norm_1 - I_norm_2) / 2^0.5;

% maximum range for mean = -1 .. 1
% maximum range for std = 0..(np/(np-1))^0.5
np = ceil(neighbors/2);
rng_m = 1.0;
rng_s = (1+1/(np-1)) ^ 0.5;

u_m = rng_m*2 / (bins-1);
bm = -rng_m:u_m:rng_m;
u_s = rng_s / (bins-1);
bs = 0:u_s:rng_s;

% f1->rng1 values in [0..bins-1]
m_high = mean(y_high,3);
m_low = mean(y_low,3);
hm_high = histc(m_high(:), bm);
hm_low = histc(m_low(:), bm);

% f2->rng2 values in [0..bins-1]
s_high = std(y_high, 0, 3);
s_low = std(y_low, 0, 3);
hs_high = histc(s_high(:), bs);
hs_low = histc(s_low(:), bs);

if strcmp(mode,'nh')
    hm_high = hm_high / sum(hm_high);
    hm_low = hm_low / sum(hm_low);
    hs_high = hs_high / sum(hs_high);
    hs_low = hs_low / sum(hs_low);
end
t=toc;
result = [hm_high hm_low hs_high hs_low];
end